05_MYB_HPV_fisher

1 Functions

library(rstatix)

2 Data

hmsc_signatue

 [1] "ANLN"         "TUBB6"        "AXL"          "IFT122"       "GFM1"         "SEPT2"        "PBRM1"        "INSIG1"       "CAPG"         "POP7"        
[11] "BHLHE40"      "MCM5"         "TRAPPC3"      "MRPS6"        "GINS2"        "NARS"         "LSM5"         "GUSB"         "EZR"          "FZD6"        
[21] "STAT2"        "TMPO"         "SLC35F5"      "MVB12A"       "FBLN1"        "COG4"         "MRPL20"       "GSDMD"        "VCL"          "MT1E"        
[31] "TUBA1C"       "ATXN10"       "USP39"        "MSH6"         "SDC2"         "OAT"          "SCAP"         "PRELID3B"     "KNTC1"        "WWTR1"       
[41] "ARF6"         "SMC5"         "DCTN2"        "AAMP"         "CSNK1A1"      "TBCB"         "KLF10"        "TPP1"         "EIF5A"        "IFT52"       
[51] "ETS2"         "EWSR1"        "NOP10"        "KRCC1"        "ASL"          "IGSF9"        "PHF5A"        "XRCC5"        "RNF168"       "UBE2T"       
[61] "IK"           "AC007250.3"   "YES1"         "SNRPD1"       "C9orf3"       "FDFT1"        "LASP1"        "SNRPC"        "PACSIN2"      "MTHFD2"      
[71] "PMVK"         "OBP2B"        "AKR1B1"       "SYPL1"        "LPAR2"        "PHB"          "THOC6"        "TMX1"         "SDHA"         "NUBP2"       
[81] "NDUFS7"       "TAOK3"        "MRPS9"        "PTRF"         "PSME2"        "MED15"        "PES1"         "RP11-19E11.1" "RP5-1102E8.2" "SEH1L"       
[91] "MCM3"         "POLR2L"       "DNAJC8"       "COPS6"        "PRKCDBP"      "KRT16"        "RNASEH2A"     "PDIA6"       

# get hpv genes and calc score
hpv_16_genes = rownames(opscc_hpvPos)[ startsWith(x =rownames(opscc_hpvPos),prefix =  "HPV16")]
hpv_score = FetchData(object = opscc_hpvPos,vars = hpv_16_genes,slot = "counts") %>% rowSums()
hpv_score_df = data.frame(hpv_score,row.names = names(hpv_score))
opscc_hpvPos = AddMetaData(object = opscc_hpvPos,metadata = hpv_score_df)

#set HPV to consider + if hpv_score>3
opscc_hpvPos$hpv = FetchData(object = opscc_hpvPos,vars = "hpv") %>%
 mutate(hpv = if_else(condition = hpv_score>=3,true = "HPV+",false = "HPV-")) %>% 
  pull(hpv)

opscc_hpvPos = ScaleData(object = opscc_hpvPos,features = hmsc_signatue)

Centering and scaling data matrix

  |                                                                                                                                                       
  |                                                                                                                                                 |   0%
  |                                                                                                                                                       
  |=================================================================================================================================================| 100%

DoHeatmap(object = opscc_hpvPos,features = hmsc_signatue,group.by = "hpv")

data = FetchData(object = opscc_hpvPos, vars = c(hmsc_signatue,"patient"))
data[] <- data %>%  group_by(patient) %>% mutate(across(is.numeric, ~ as.numeric(scale(.)))) %>% ungroup() 
data = data %>% select(-patient) 
data = data %>% t() %>% as.data.frame()

annotation_data = FetchData(object = opscc_hpvPos, vars = c("patient","hpv")) %>% dplyr::arrange(patient,hpv)
column_ha = HeatmapAnnotation(hpv = annotation_data[,2])
data = data[,rownames(annotation_data)] #order data like annotation
data = na.omit(data)
p = ComplexHeatmap::Heatmap(
  data,
  show_column_names = F,
  row_names_gp = grid::gpar(fontsize = 7),
  cluster_rows = T,
  top_annotation = column_ha,
  name = "Z-score expression",use_raster = F,cluster_columns = F,column_split   = annotation_data$patient,
  column_gap = unit(2, "mm"), border = TRUE,show_parent_dend_line = FALSE, show_column_dend = FALSE,cluster_column_slices = F
)

Warning: The input is a data frame, convert it to a matrix.The automatically generated colors map from the minus and plus 99^th of the absolute values in the matrix. There are outliers in the matrix
whose patterns might be hidden by this color mapping. You can manually set the color to `col` argument.

Use `suppressMessages()` to turn off this message.

  
p

# add score
hpv_signature_score = FetchData(object = opscc_hpvPos,vars = hmsc_signatue,slot = "data") %>% rowMeans()

Warning: The following requested variables were not found: AC007250.3, RP5-1102E8.2

hpv_signature = data.frame(hpv_signature_score,row.names = names(hpv_signature_score))
opscc_hpvPos = AddMetaData(object = opscc_hpvPos,metadata = hpv_signature)

genes_by_tp = FetchData(object = opscc_hpvPos,vars = c("hpv", "hpv_signature_score","patient"),slot = "data") 



#plot and split by patient:
stat.test <- genes_by_tp %>%
  group_by(patient) %>% 
  # filter (patient != "OP13") %>% 
  t_test(hpv_signature_score ~ hpv) %>%
  adjust_pvalue(method = "fdr") %>%
  add_significance() %>%
  add_xy_position(x = "hpv") %>%
  mutate(p.adj = signif(p.adj, digits = 2))
stat.test

plt = ggplot(genes_by_tp, aes(x = hpv, y = hpv_signature_score)) +
  geom_violin(scale = "width",aes(fill = hpv)) +
  geom_boxplot(width = 0.2, outlier.shape = NA) +
  stat_pvalue_manual(stat.test, label = "{p.adj}")  + #add p value
  facet_wrap("patient")+
  scale_y_continuous(expand = expansion(mult = c(0.05, 0.2))) + # Add 10% spaces between the p-value labels and the plot border
  ggtitle("HMSC HPV signature")
plt